Library preperation

library(dplyr)
library(GenomicFeatures)
library(readr)
library(enrichR)
library(DESeq2)
library(ensembldb)
library(biomaRt)
library(clusterProfiler)
library(DOSE)

library(RColorBrewer)
library(gplots)
library(plotly)
library(pheatmap)
library(ggplot2)
library(easylabel)
library(enrichplot)
library(ggrepel)

library(AnnotationDbi)
library(org.Hs.eg.db)
library(pathview)
library(gage)
library(gageData)
library(tximport)
library(txdbmaker)

Get nnotation data

txdb <- makeTxDbFromGFF("gencode.v46.annotation.gff3")
## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ...
## Warning in .get_cds_IDX(mcols0$type, mcols0$phase): The "phase" metadata column contains non-NA values for features of type
##   stop_codon. This information was ignored.
## OK
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(txdb, k, "GENEID", "TXNAME")
## 'select()' returned 1:1 mapping between keys and columns
head(tx2gene)

Experiment summary & quantification

samples <- read.table("samples.txt",sep='\t', header = TRUE)
to_drop <- c("BIOMATERIAL_PROVIDER","BIOMATERIAL_PROVIDER","Center.Name","DATASTORE.provider", "Consent", "Bytes", "DATASTORE.filetype", "DATASTORE.region", "LibrarySelection", "Collection_Date", "geo_loc_name_country", "geo_loc_name_country_continent", "geo_loc_name","Assay.Type","BioSampleModel","Instrument", "Platform", "Organism", "version", "SRA.Study")
columns_to_keep <- setdiff(names(samples), to_drop)
samples <- samples %>% dplyr::select(all_of(columns_to_keep)) |> filter(LibrarySource =="TRANSCRIPTOMIC")

samples$tissue <- gsub("tibialis anterior \\(TA\\)", "TA", samples$tissue)
samples$tissue <- gsub("vastus lateralis \\(VL\\)", "VL", samples$tissue)

samples$tissue <- as.factor(samples$tissue)
samples$health_state <- as.factor(samples$health_state)
samples$sex <- as.factor(samples$sex)
samples$BulkRNAseq_sequencing_batch <- as.factor(samples$BulkRNAseq_sequencing_batch)
summary_df <- samples %>% filter(health_state=='healthy') %>%
  group_by(BulkRNAseq_sequencing_batch, sex, tissue) %>%
  summarize(Sample_Count = n())
## `summarise()` has grouped output by 'BulkRNAseq_sequencing_batch', 'sex'. You
## can override using the `.groups` argument.
samples_TA <- samples |> dplyr::filter(tissue == 'TA') #|> filter(sex=="male")
files <- file.path("quants", samples_TA$Run, "quant.sf")
names(files) <- paste0(samples_TA$Run)
txi.ta <- tximport(files, type = "salmon", tx2gene = tx2gene)
## reading in files with read_tsv
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 
## summarizing abundance
## summarizing counts
## summarizing length
head(txi.ta$counts,3)
##                    SRR24750870 SRR24750872 SRR24750875 SRR24750884 SRR24750886
## ENSG00000000003.16      37.000      39.000     120.190         104      71.000
## ENSG00000000005.6       26.000       1.000       8.000           0       4.000
## ENSG00000000419.14     277.001     308.001     902.406         465     572.387
##                    SRR24750888 SRR24750890 SRR24750891 SRR24750893 SRR24750895
## ENSG00000000003.16     102.042      59.000       64.00      80.999      66.000
## ENSG00000000005.6        4.000       4.000        1.00       5.000       1.000
## ENSG00000000419.14     666.009     351.855      567.08     316.001     322.999
##                    SRR24750897 SRR24750899 SRR24750902 SRR25143906 SRR25143907
## ENSG00000000003.16      97.000      85.001      37.000     167.001     144.098
## ENSG00000000005.6       11.000       4.000      46.000     247.415     247.762
## ENSG00000000419.14     302.001     422.999     284.343     722.819     584.200
##                    SRR25143908 SRR25143909 SRR25143910 SRR25143911 SRR25143912
## ENSG00000000003.16     127.190     178.001      65.999     151.054     114.000
## ENSG00000000005.6      409.708      12.000       5.000      16.000      47.000
## ENSG00000000419.14     619.122     730.491     431.999     778.456     553.069
##                    SRR25143913
## ENSG00000000003.16     133.000
## ENSG00000000005.6        8.000
## ENSG00000000419.14     614.075

Run DESeq2

ddsTxi <- DESeqDataSetFromTximport(txi.ta,
                                   colData = samples_TA,
                                   design = ~ health_state + BulkRNAseq_sequencing_batch + sex )
## using counts and average transcript lengths from tximport
keep <- rowSums(counts(ddsTxi)) >= 10
ddsTxi <- ddsTxi[keep,]
dds <- DESeq(ddsTxi)
## estimating size factors
## using 'avgTxLength' from assays(dds), correcting for library size
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 150 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
res <- results(dds,contrast=c("health_state","DMD","healthy"))
saveRDS(dds, 'dds_ta.RDS')
summary(res)
## 
## out of 29429 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 5091, 17%
## LFC < 0 (down)     : 4577, 16%
## outliers [1]       : 91, 0.31%
## low counts [2]     : 3986, 14%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
res <- res[order(-abs(res$log2FoldChange)), ]
resDF <- as.data.frame(res)
resDF <- dplyr::as_tibble(resDF, rownames = "ensembl_gene_id")

resFiltered <- resDF %>% dplyr::filter(abs(log2FoldChange)>1, padj<0.05)
resFiltered <- dplyr::as_tibble(resFiltered, rownames = "ensembl_gene_id")

Graphical inspection

plotDispEsts(dds, main="Dispersion plot")

rld <- rlogTransformation(dds)
head(assay(rld))
##                    SRR24750870 SRR24750872 SRR24750875 SRR24750884 SRR24750886
## ENSG00000000003.16    6.329671    6.051057    6.267613    6.779514    6.176066
## ENSG00000000005.6     4.559225    2.904288    3.349846    2.762746    3.672551
## ENSG00000000419.14    8.591156    8.770587    9.307182    8.997897    9.203417
## ENSG00000000457.14    7.589054    7.869972    8.058020    8.274363    7.901477
## ENSG00000000460.17    6.126614    5.272208    5.909147    6.000162    5.593122
## ENSG00000000938.13    5.477545    6.233680    6.394446    5.567846    8.023214
##                    SRR24750888 SRR24750890 SRR24750891 SRR24750893 SRR24750895
## ENSG00000000003.16    6.244494    6.493658    6.119276    6.762724    6.592733
## ENSG00000000005.6     3.090547    3.234220    2.902277    3.313863    2.884141
## ENSG00000000419.14    9.171032    8.740413    9.324809    8.656339    8.704264
## ENSG00000000457.14    8.053547    7.691997    7.808946    7.911316    7.842755
## ENSG00000000460.17    6.169462    5.536892    6.231253    5.722466    6.324115
## ENSG00000000938.13    5.491031    5.875441    5.238462    5.942633    5.990384
##                    SRR24750897 SRR24750899 SRR24750902 SRR25143906 SRR25143907
## ENSG00000000003.16    6.489076    6.625695    6.098108    6.397513    6.927003
## ENSG00000000005.6     3.969479    3.255026    5.183687    6.139530    6.525705
## ENSG00000000419.14    8.608024    8.932316    8.694287    8.795892    8.952481
## ENSG00000000457.14    7.737008    7.944524    8.048814    7.972602    7.690648
## ENSG00000000460.17    6.052215    5.322183    5.690707    6.500251    6.008938
## ENSG00000000938.13    6.434979    6.320757    5.837344    6.774240    6.796708
##                    SRR25143908 SRR25143909 SRR25143910 SRR25143911 SRR25143912
## ENSG00000000003.16    6.451050    6.820612    6.406500    6.557796    6.592475
## ENSG00000000005.6     6.840662    3.447859    3.260614    3.625384    4.821752
## ENSG00000000419.14    8.936241    8.954982    9.018376    9.017794    8.982253
## ENSG00000000457.14    7.744739    7.939098    7.941927    7.745713    7.639098
## ENSG00000000460.17    6.421365    6.514294    6.369991    6.308550    6.132154
## ENSG00000000938.13    6.304516    6.216860    6.942990    6.475628    5.669449
##                    SRR25143913
## ENSG00000000003.16    6.602341
## ENSG00000000005.6     3.317934
## ENSG00000000419.14    8.879044
## ENSG00000000457.14    8.056516
## ENSG00000000460.17    6.420569
## ENSG00000000938.13    6.416174
mycols <- brewer.pal(8, "Dark2")[1:length(unique(samples_TA$health_state))]
sampleDists <- as.matrix(dist(t(assay(rld))))
heatmap.2(as.matrix(sampleDists), key=F, trace="none",
          col=colorpanel(100, "black", "white"),
          ColSideColors=mycols[samples_TA$health_state],
          RowSideColors=mycols[samples_TA$health_state],
          margin=c(10, 10), main="Sample Distance Matrix")
legend(x=0,y=1.1, legend = levels(samples_TA$health_state), fill = mycols, title = "Categories",bty='n')

selected <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("health_state", "sex")])
pheatmap(assay(rld)[selected,], cluster_rows=FALSE, show_rownames=TRUE,
         cluster_cols=FALSE, annotation_col=df)

DESeq2::plotPCA(rld, intgroup=c("health_state"))
## using ntop=500 top features by variance

pcaData <- plotPCA(rld, intgroup=c("health_state",'sex'), returnData=TRUE)
## using ntop=500 top features by variance
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=health_state, shape=sex)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()

# hist(res$pvalue, breaks=50, col="grey")
# DESeq2::plotMA(dds, ylim=c(-1,1))
# 
# # Volcano plot
# with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2)))
# with(subset(res, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red"))
fig <- plot_ly(x = res$pvalue, type = "histogram")
fig <- fig %>% layout(title = "Histogram of p-values",
                     xaxis = list(title = "Value"),
                     yaxis = list(title = "Frequency"))
fig
## Warning: Ignoring 91 observations
easyMAplot(res, output_shiny=FALSE)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
res <- res[order(-abs(res$log2FoldChange)), ]
resFiltered <- resDF %>% dplyr::filter(abs(log2FoldChange)>1.5,padj<0.05)

resDF$ensembl_gene_id <- gsub("\\.\\d+$", "", resDF$ensembl_gene_id)
gene_info <- getBM(c("ensembl_gene_id","hgnc_symbol","entrezgene_id","description","entrezgene_description"), "ensembl_gene_id", resDF$ensembl_gene_id, ensembl) |> group_by(ensembl_gene_id) |> dplyr::slice_head(n = 1) |>  ungroup()
resDF <- resDF %>%
  left_join(gene_info, by = "ensembl_gene_id")
res1111 <- resDF %>% dplyr::filter(padj<0.05)
top_peaks <- as.data.frame(head(resDF,10))
a <- list()
for (i in seq_len(nrow(top_peaks))) {
  m <- top_peaks[i, ]
  annot = ifelse(m[["hgnc_symbol"]]=="", "novel", m[["hgnc_symbol"]])
  ax <- sample(c(-1, 1), 1) * runif(1, 5, 60)
  ay <- sample(c(-1, 1), 1) * runif(1, 5, 40)
  a[[i]] <- list(
    x = m[["log2FoldChange"]],
    y = -log10(m[["pvalue"]]),
    text = annot,
    xref = "x",
    yref = "y",
    showarrow = TRUE,
    arrowhead = 0.5,
    ax = ax,
    ay = ay
  )
}
easyVolcano(res, fccut = 2, output_shiny=FALSE) %>%   layout(annotations = a)
## Warning: Ignoring 91 observations
data(kegg.sets.hs)
data(sigmet.idx.hs)
kegg.sets.hs <- kegg.sets.hs[sigmet.idx.hs]
foldchanges <- resDF$log2FoldChange
names(foldchanges) <- resDF$entrezgene_id
keggres <- gage(foldchanges, gsets=kegg.sets.hs, same.dir=FALSE)
lapply(keggres, head)
## $greater
##                                                 p.geomean stat.mean
## hsa04512 ECM-receptor interaction            0.0003702404  3.441015
## hsa04610 Complement and coagulation cascades 0.0005943481  3.328515
## hsa04974 Protein digestion and absorption    0.0017050649  2.979736
## hsa04514 Cell adhesion molecules (CAMs)      0.0021366136  2.884320
## hsa04115 p53 signaling pathway               0.0403289882  1.760560
## hsa04145 Phagosome                           0.0629244296  1.535402
##                                                     p.val      q.val set.size
## hsa04512 ECM-receptor interaction            0.0003702404 0.04754784       85
## hsa04610 Complement and coagulation cascades 0.0005943481 0.04754784       56
## hsa04974 Protein digestion and absorption    0.0017050649 0.08546454       70
## hsa04514 Cell adhesion molecules (CAMs)      0.0021366136 0.08546454      123
## hsa04115 p53 signaling pathway               0.0403289882 1.00000000       67
## hsa04145 Phagosome                           0.0629244296 1.00000000      138
##                                                      exp1
## hsa04512 ECM-receptor interaction            0.0003702404
## hsa04610 Complement and coagulation cascades 0.0005943481
## hsa04974 Protein digestion and absorption    0.0017050649
## hsa04514 Cell adhesion molecules (CAMs)      0.0021366136
## hsa04115 p53 signaling pathway               0.0403289882
## hsa04145 Phagosome                           0.0629244296
## 
## $stats
##                                              stat.mean     exp1
## hsa04512 ECM-receptor interaction             3.441015 3.441015
## hsa04610 Complement and coagulation cascades  3.328515 3.328515
## hsa04974 Protein digestion and absorption     2.979736 2.979736
## hsa04514 Cell adhesion molecules (CAMs)       2.884320 2.884320
## hsa04115 p53 signaling pathway                1.760560 1.760560
## hsa04145 Phagosome                            1.535402 1.535402
keggrespathwaysUp <- data.frame(id=rownames(keggres$greater), keggres$greater) %>%
  tibble::as_tibble() %>%
  dplyr::filter(row_number()<=6) %>%
  .$id %>%
  as.character()
keggrespathwaysUp
## [1] "hsa04512 ECM-receptor interaction"           
## [2] "hsa04610 Complement and coagulation cascades"
## [3] "hsa04974 Protein digestion and absorption"   
## [4] "hsa04514 Cell adhesion molecules (CAMs)"     
## [5] "hsa04115 p53 signaling pathway"              
## [6] "hsa04145 Phagosome"
keggresids <- substr(c(keggrespathwaysUp), start=1, stop=8)
keggresids
## [1] "hsa04512" "hsa04610" "hsa04974" "hsa04514" "hsa04115" "hsa04145"
plot_pathway <- function(pid) pathview(gene.data=foldchanges, pathway.id=pid, species="hsa", new.signature=FALSE)
#detach("package:dplyr", unload=T)
tmp <- sapply(keggresids, function(pid) pathview(gene.data=foldchanges, pathway.id=pid, species="hsa"))
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory D:/VLTA
## Info: Writing image file hsa04512.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory D:/VLTA
## Info: Writing image file hsa04610.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory D:/VLTA
## Info: Writing image file hsa04974.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory D:/VLTA
## Info: Writing image file hsa04514.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory D:/VLTA
## Info: Writing image file hsa04115.pathview.png
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory D:/VLTA
## Info: Writing image file hsa04145.pathview.png
resDF <- resDF %>% mutate(diffexpressed = dplyr::case_when(
  log2FoldChange > 0 & padj < 0.05 ~ 'UP',
  log2FoldChange < 0 & padj < 0.05 ~ 'DOWN',
  padj > 0.05 ~ 'NO'
))
original_gene_list <- resDF$log2FoldChange
names(original_gene_list) <- resDF$ensembl_gene_id
gene_list<-na.omit(original_gene_list)
gene_list = sort(gene_list, decreasing = TRUE)
gse <- gseGO(geneList=gene_list, 
             ont ="ALL", 
             keyType = "ENSEMBL", 
             minGSSize = 3, 
             maxGSSize = 800, 
             pvalueCutoff = 0.05, 
             verbose = TRUE, 
             OrgDb = org.Hs.eg.db, 
             pAdjustMethod = "none")
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (12.01% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...
dotplot(gse, showCategory=10, font.size=7, split=".sign") + facet_grid(.~.sign)

#emapplot(gse, showCategory = 10)
websiteLive <- getOption("enrichR.live")
if (websiteLive) {
    listEnrichrSites()
    setEnrichrSite("Enrichr") # Human genes   
}
## Enrichr ... Connection is Live!
## FlyEnrichr ... Connection is Live!
## WormEnrichr ... Connection is Live!
## YeastEnrichr ... Connection is Live!
## FishEnrichr ... Connection is Live!
## OxEnrichr ... Connection is Live!
## Connection changed to https://maayanlab.cloud/Enrichr/
## Connection is Live!
if (websiteLive) dbs <- listEnrichrDbs()
dbs <- c("GO_Molecular_Function_2023", "GO_Cellular_Component_2023", "GO_Biological_Process_2023","GO_Molecular_Function_2018", "GO_Cellular_Component_2018", "GO_Biological_Process_2018")
if (websiteLive) {
    enriched <- enrichr(names(gene_list), dbs)
}
## Uploading data to Enrichr... Done.
##   Querying GO_Molecular_Function_2023... Done.
##   Querying GO_Cellular_Component_2023... Done.
##   Querying GO_Biological_Process_2023... Done.
##   Querying GO_Molecular_Function_2018... Done.
##   Querying GO_Cellular_Component_2018... Done.
##   Querying GO_Biological_Process_2018... Done.
## Parsing results... Done.
df <- resDF |> dplyr::filter(padj<0.05)
up_regulated <- df |> dplyr::filter(log2FoldChange > 0) |>   dplyr::mutate(ensembl_gene_id = gsub("\\.\\d+$", "", ensembl_gene_id)) |> dplyr::select(ensembl_gene_id)
write.table(up_regulated, 'up_DMD_healthy.txt', row.names = FALSE, col.names = FALSE, quote = FALSE)
down_regulated <- df |> dplyr::filter(log2FoldChange < 0) |>   dplyr::mutate(ensembl_gene_id = gsub("\\.\\d+$", "", ensembl_gene_id)) |> dplyr::select(ensembl_gene_id)
write.table(down_regulated, 'down_DMD_healthy.txt', row.names = FALSE, col.names = FALSE, quote = FALSE)
# Choose databases
dbs <- c("GO_Biological_Process_2021", "KEGG_2021_Human", "WikiPathways_2019_Human", "Jensen_DISEASES")
genes <- names(gene_list)
# Run enrichment analysis 
enriched_results <- enrichr(genes, dbs)
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
##   Querying KEGG_2021_Human... Done.
##   Querying WikiPathways_2019_Human... Done.
##   Querying Jensen_DISEASES... Done.
## Parsing results... Done.
# Access results for a specific database
go_results <- enriched_results[["GO_Biological_Process_2021"]]
kegg_results <- enriched_results[["KEGG_2021_Human"]]
egoCC <- enrichGO(gene          = as.character(gene_info$entrezgene_id),
                #universe      = as.character(universe$entrezgene_id),
                OrgDb         = org.Hs.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)
egoCCDF <- egoCC@result
write.table(egoCCDF, 'CC_DMD_H.csv')
egoBP <- enrichGO(gene          = as.character(gene_info$entrezgene_id),
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)
egoBPBP <- egoBP@result
write.table(egoBPBP, 'BP_DMD_H.csv')
sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Polish_Poland.utf8  LC_CTYPE=Polish_Poland.utf8   
## [3] LC_MONETARY=Polish_Poland.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Polish_Poland.utf8    
## 
## time zone: Europe/Warsaw
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] txdbmaker_1.0.0             tximport_1.32.0            
##  [3] gageData_2.42.0             gage_2.54.0                
##  [5] pathview_1.44.0             org.Hs.eg.db_3.19.1        
##  [7] ggrepel_0.9.5               enrichplot_1.24.0          
##  [9] easylabel_0.2.8             pheatmap_1.0.12            
## [11] plotly_4.10.4               ggplot2_3.5.1              
## [13] gplots_3.1.3.1              RColorBrewer_1.1-3         
## [15] DOSE_3.30.1                 clusterProfiler_4.12.0     
## [17] biomaRt_2.60.0              ensembldb_2.28.0           
## [19] AnnotationFilter_1.28.0     DESeq2_1.44.0              
## [21] SummarizedExperiment_1.34.0 MatrixGenerics_1.16.0      
## [23] matrixStats_1.3.0           enrichR_3.2                
## [25] readr_2.1.5                 GenomicFeatures_1.56.0     
## [27] AnnotationDbi_1.66.0        Biobase_2.64.0             
## [29] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [31] IRanges_2.38.0              S4Vectors_0.42.0           
## [33] BiocGenerics_0.50.0         dplyr_1.1.4                
## 
## loaded via a namespace (and not attached):
##   [1] later_1.3.2              splines_4.4.1            BiocIO_1.14.0           
##   [4] bitops_1.0-7             ggplotify_0.1.2          filelock_1.0.3          
##   [7] tibble_3.2.1             polyclip_1.10-6          graph_1.82.0            
##  [10] XML_3.99-0.16.1          lifecycle_1.0.4          httr2_1.0.1             
##  [13] lattice_0.22-6           MASS_7.3-60.2            crosstalk_1.2.1         
##  [16] magrittr_2.0.3           sass_0.4.9               rmarkdown_2.27          
##  [19] jquerylib_0.1.4          yaml_2.3.8               httpuv_1.6.15           
##  [22] cowplot_1.1.3            DBI_1.2.3                abind_1.4-5             
##  [25] zlibbioc_1.50.0          purrr_1.0.2              ggraph_2.2.1            
##  [28] RCurl_1.98-1.14          yulab.utils_0.1.4        WriteXLS_6.6.0          
##  [31] tweenr_2.0.3             rappdirs_0.3.3           GenomeInfoDbData_1.2.12 
##  [34] tidytree_0.4.6           codetools_0.2-20         DelayedArray_0.30.1     
##  [37] DT_0.33                  xml2_1.3.6               ggforce_0.4.2           
##  [40] tidyselect_1.2.1         aplot_0.2.3              UCSC.utils_1.0.0        
##  [43] farver_2.1.2             viridis_0.6.5            BiocFileCache_2.12.0    
##  [46] GenomicAlignments_1.40.0 jsonlite_1.8.8           tidygraph_1.3.1         
##  [49] tools_4.4.1              progress_1.2.3           treeio_1.28.0           
##  [52] Rcpp_1.0.12              glue_1.7.0               gridExtra_2.3           
##  [55] SparseArray_1.4.8        xfun_0.44                qvalue_2.36.0           
##  [58] withr_3.0.0              fastmap_1.2.0            fansi_1.0.6             
##  [61] caTools_1.18.2           digest_0.6.35            mime_0.12               
##  [64] R6_2.5.1                 gridGraphics_0.5-1       colorspace_2.1-0        
##  [67] GO.db_3.19.1             gtools_3.9.5             RSQLite_2.3.7           
##  [70] utf8_1.2.4               tidyr_1.3.1              generics_0.1.3          
##  [73] data.table_1.15.4        rtracklayer_1.64.0       htmlwidgets_1.6.4       
##  [76] prettyunits_1.2.0        graphlayouts_1.1.1       httr_1.4.7              
##  [79] S4Arrays_1.4.1           scatterpie_0.2.3         pkgconfig_2.0.3         
##  [82] gtable_0.3.5             blob_1.2.4               XVector_0.44.0          
##  [85] shadowtext_0.1.3         htmltools_0.5.8.1        fgsea_1.30.0            
##  [88] ProtGenerics_1.36.0      scales_1.3.0             png_0.1-8               
##  [91] ggfun_0.1.5              knitr_1.47               rstudioapi_0.16.0       
##  [94] tzdb_0.4.0               reshape2_1.4.4           rjson_0.2.21            
##  [97] nlme_3.1-164             curl_5.2.1               cachem_1.1.0            
## [100] stringr_1.5.1            KernSmooth_2.23-24       shinycssloaders_1.0.0   
## [103] parallel_4.4.1           HDO.db_0.99.1            restfulr_0.0.15         
## [106] pillar_1.9.0             grid_4.4.1               vctrs_0.6.5             
## [109] promises_1.3.0           dbplyr_2.5.0             xtable_1.8-4            
## [112] Rgraphviz_2.48.0         KEGGgraph_1.64.0         evaluate_0.24.0         
## [115] cli_3.6.2                locfit_1.5-9.9           compiler_4.4.1          
## [118] Rsamtools_2.20.0         rlang_1.1.4              crayon_1.5.3            
## [121] labeling_0.4.3           plyr_1.8.9               fs_1.6.4                
## [124] stringi_1.8.4            viridisLite_0.4.2        BiocParallel_1.38.0     
## [127] munsell_0.5.1            Biostrings_2.72.1        lazyeval_0.2.2          
## [130] GOSemSim_2.30.0          Matrix_1.7-0             hms_1.1.3               
## [133] patchwork_1.2.0          bit64_4.0.5              shiny_1.8.1.1           
## [136] KEGGREST_1.44.1          highr_0.11               igraph_2.0.3            
## [139] memoise_2.0.1            bslib_0.7.0              ggtree_3.12.0           
## [142] fastmatch_1.1-4          bit_4.0.5                shinybusy_0.3.3         
## [145] ape_5.8                  gson_0.1.0